%% "Econometrics of the Hodrick-Prescott Filter," Robert de Jong and Neslihan Sakarya (2015)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% This file allows you to compare
 %(1) The numerical equivalence of the function f_T_lambda (given in Theorem 1 of the paper) and its 3rd order difference formulation (given in Lemma 4 of the paper)

 %(2) The numerical equivalence of the function g_T_lambda (given in Theorem 1 of the paper) and its 3rd order difference formulation (given in Lemma 4 of the paper)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clear all;
close all;
clc;

%% The equivalence of f and g functions with their 3 order difference formulation are numerically shown
%First, we establish the equivalence for f and then g functions

T=10; %the total number of observations
lm=1; %lambda -the smoothing parameter

%% ==============================================================%
%Showing the Equivalence of f function with its 3rd Order Difference Formula
%==============================================================%


%% Original f

f=zeros(T,1);

for j=2:T;
    for m=1:T;
        
        f(m)=f(m)+(1/T)*cos(pi*m*(j-1)/T)/(1+16*lm*sin(pi*(j-1)/(2*T))^4);
    end
end
for m=1:T;
    f(m)=f(m)+(1/(2*T))+cos(pi*m)*(1/(2*T))*(1/(1+16*lm));
end

%% Third Order Difference formula of f
% First, define h(j/T)
h=zeros(T,1);
for j=1:T-1;
    
    h(j)=1/(1+16*lm*sin(pi*j/(2*T))^4);
end

% Secondly, calculate the last component of the formula 
% i.e. sum j=2 to T-3 D^3h(.)*sin(.) 
ss=zeros(T,1);

for m=1:T;
    for j=2:T-3;
    
    ss(m)=ss(m)+(h(j+2)-3*h(j+1)+3*h(j)-h(j-1))*sin(pi*m*(j+(1/2))/T);
end
end

%Writing the third order difference formula explicitly for different values
%of m

f3=zeros(T,1);
for m=1:T;
    
    f3(m)=(1/(2 *T))+cos(pi*m)*(1/(2*T))*(1/(1+16*lm))-(1/T)*(1/(2*sin(pi*m/(2*T))))*h(1)*sin(pi*m/(2*T))+(1/T)*(1/(2*sin(pi*m/(2*T))))*h(T-1)*sin(pi*m*(T-(1/2))/T)-(1/T)*(1/(2*sin(pi*m/(2*T))))^2*(h(2)-h(1))*cos(pi*m/T)+(1/T)*(1/(2*sin(pi*m/(2*T))))^2*(h(T-1)-h(T-2))*cos(pi*m*(T-1)/T)+(1/T)*(1/(2*sin(pi*m/(2*T))))^3*(h(3)-2*h(2)+h(1))*sin(pi*m*3/(2*T))-(1/T)*(1/(2*sin(pi*m/(2*T))))^3*(h(T-1)-2*h(T-2)+h(T-3))*sin(pi*m*(T-(3/2))/T)+(1/T)*(1/(2*sin(pi*m/(2*T))))^3*ss(m);
end


%% ==============================================================
%%%%%%%%%Showing the Equivalence of g function with its 3rd Order
%%%%%%%%%Difference Formula
%==============================================================

%% Original g function

ht=zeros(T-1,1); %% Construction of h_tilda

for j=1:T-1;
    
    ht(j)=(sin(pi*j/(2*T))^2)/(1+16*lm*(sin(pi*j/(2*T)))^4);
end

M1=zeros(T,T-1);
M2=zeros(T-1,1);

for m=2:T;
    for j=1:T-1;
        
        M1(m,j)=cos(pi*j*(m-(1/2))/T);
        M2(j)=cos(pi*j/(2*T))*ht(j);
    end
end

g=sqrt(2)*M1*M2/T;

%% Third order difference of g 

% To use short cut I define k=the sin func's in the sum  
k=zeros(T,T-3);

for j=2:T-3; %%% j=1 ile basliyordu
   for m=2:T;
    
       
    k(m,j)=((sin(pi*m*(j+(1/2))/T))/(2*sin(pi*m/(2*T)))^3)+((sin(pi*(m-1)*(j+(1/2))/T))/(2*sin(pi*(m-1)/(2*T)))^3);
    

   end
end

% Third difference equations

D3h=zeros(T-3,1);
D3h(1)=ht(3)-(3*ht(2))+(3*ht(1)); %1st row of D3h
for j=2:T-3;
    
    D3h(j)=ht(j+2)-(3*ht(j+1))+(3*ht(j))-ht(j-1);
end

% S is the last line of the g function ---sum from j=2, T-3
S=k*D3h;

% The Third Order Difference formula of g function

c=sqrt(2)/(2*T);

gn=zeros(T,1);
%Here we do not calculate gn(1) because the function is undefined
%(i.e. one of the denominator term sin(pi(m-1)/(2T))=0) when m=1
for m=2:T;
    
   gn(m)=-c*ht(1)-c*(ht(2)-ht(1))*((cos(pi*m/T)/(2*sin(pi*m/(2*T)))^2)+(cos(pi*(m-1)/T)/(2*sin(pi*(m-1)/(2*T)))^2))+c*(ht(T-1)-ht(T-2))*((cos(pi*m*(T-1)/T)/(2*sin(pi*m/(2*T)))^2)+(cos(pi*(m-1)*(T-1)/T)/(2*sin(pi*(m-1)/(2*T)))^2))+c*(ht(3)-(2*ht(2))+ht(1))*((sin(3*pi*m/(2*T))/(2*sin(pi*m/(2*T)))^3)+(sin(3*pi*(m-1)/(2*T))/(2*sin(pi*(m-1)/(2*T)))^3))-c*(ht(T-1)-(2*ht(T-2))+ht(T-3))*((sin(pi*m*(T-(3/2))/T)/(2*sin(pi*m/(2*T)))^3)+(sin(pi*(m-1)*(T-(3/2))/T)/(2*sin(pi*(m-1)/(2*T)))^3))+c*S(m);

end

%% The results below will be displayed on the command window

% The result for f

Diff=f-f3; %The difference between the orginal f and its 3rd order difference formula
Rf=[f, f3, Diff];
disp('    f_T_lambda        3rd Order Dif Form.   Difference')
disp(Rf)

% The result for g

Difg=gn-g; %The difference between the orginal g and its 3rd order difference formula
Rg=[g, gn, Difg];
disp('   g_T_lambda        3rd Order Dif Form.    Difference')
disp(Rg)


    
    